#include "2dpet.h"

static float image[196 * 128], image1[196 * 128];

// Print TGA header
void PrintTGAHeader( FILE * ofile, int width, int height ) {
	fputc(0,ofile);
	fputc(0,ofile);
	fputc(2,ofile);
	for(int i=3;i<12;i++) fputc(0,ofile);
	fputc(width % 256, ofile);
	fputc(width / 256, ofile);
	fputc(height % 256, ofile);
	fputc(height / 256, ofile);
	fputc(24,ofile);
	fputc(32,ofile);
}

// Print TGA header
void PrintTGAColor( FILE * ofile, double value, double gamma = 0.8 ) {
	int r, g, b;
	if (value < 0) value = 0;
	if (value > 1) value = 1;
	r = g = b = (int)((1 - pow(value, gamma)) * 255.5);
	fputc(b, ofile); fputc(g, ofile); fputc(r, ofile);
}

void Message(char * buffer, bool up) {    
	glColor3f(1, 1, 1);
	glRasterPos2f( 0.01, (up) ? 0.85 : -0.98 );
	for( int i = 0; i < 80; i++ ) {
		if ( buffer[i] == '\0' ) break;
		else glutBitmapCharacter( GLUT_BITMAP_HELVETICA_10, buffer[i] );
	}
}

void Voxels :: PrintImage( char * ofname, double maxvoxel, bool bar ) {
	cout << "Image: " << ofname << '\n';
	static char buffer[80];
	sprintf(buffer,"%s.ima", ofname);
	PrintPlot( buffer, maxvoxel );

	sprintf(buffer,"%s.tga", ofname);
	FILE * ofile = fopen(buffer, "wb");
	int Xres = (bar) ? xres : yres;
	PrintTGAHeader( ofile, Xres, yres );

	for ( int i = 0; i < yres; ++i ) {
		for ( unsigned int j = 0; j < yres; ++j ) {
			image[i * vrows * 4 + j] = at(i / image_scale, j / image_scale) / maxvoxel;
		}
	}
	glClear(GL_COLOR_BUFFER_BIT);		
	glRasterPos2f(-1, -1);
	glDrawPixels(yres, yres, GL_LUMINANCE, GL_FLOAT, image); 		
	sprintf(buffer,"%d", 0);    Message(buffer, false);
	sprintf(buffer,"%d", (int)maxvoxel); Message(buffer, true);
	glRasterPos2f(-1, -1);
	glReadPixels(0, 0, Xres, yres, GL_LUMINANCE, GL_FLOAT, image); 

	for ( int i =  yres - 1; i >= 0; --i ) {
		for ( unsigned int j = 0; j < Xres; ++j ) {
			if (bar && j >= yres - 3 && j < yres) image[i * Xres + j] = (double)(i - 1)/yres;
			PrintTGAColor( ofile, (double)image[i * Xres + j], 1 );
		}
	}
	fclose(ofile);
}

void Voxels :: PrintPlot( char * ofname, double maxvoxel ) {
	return;

	cout << "\n Image: " << ofname << '\n';
	FILE * ofile = fopen(ofname, "w");

	for ( int i = 0; i < vrows; ++i ) {
		for ( unsigned int j = 0; j < vrows; ++j ) {
			if (i == vrows/4 && j == vrows/4) fprintf(ofile, "%f ", maxvoxel);
			else						      fprintf(ofile, "%f ", at(i, j) < maxvoxel ? at(i, j) : maxvoxel);
		}
		fprintf(ofile, "\n");
	}
	fclose(ofile);
}

void Lors :: PrintSinogram( char * ofname, double maxvalue, bool bar ) {
	cout << "\n Sinogram: " << ofname << '\n';
	static char buffer[80];
	sprintf(buffer,"%s.tga", ofname);
	FILE * tgafile = fopen(buffer, "wb");
	sprintf(buffer,"%s.ima", ofname);
	FILE * imafile = fopen(buffer, "w");

	int Xres = (bar) ? xres : yres;
	PrintTGAHeader( tgafile, Xres, yres );

	for(int i = 0; i < ndets; i++ ) {
		for(int j = 0; j < ndets; j++) {
			double alpha = M_PI * (i + 0.5) / ndets;
			double displace = (((double)j + 0.5) / ndets - 0.5) * volume;
			double dphi = 2.0 * acos(fabs(displace) / R);
			double phi = alpha - dphi / 2.0;
			if (displace < 0) phi += M_PI;
			double phi1 = phi + dphi;
			if (phi < 0) phi += 2 * M_PI;
			if (phi1 >= 2 * M_PI) phi1 -= 2 * M_PI;
			if (phi > phi1) {
				double tmp = phi;
				phi = phi1;
				phi1 = tmp;
			}

			int iphi = (int)(ndets * phi / 2.0 / M_PI);
			if (iphi > ndets-1) iphi = ndets - 1; 
			int idphi = (int)(ndets * (phi1 - phi) / 2.0 / M_PI + 0.5);
			if (iphi + idphi > ndets-1) idphi = ndets + 1 - iphi;
			double val = at(iphi, idphi) / maxvalue;
			if (idphi <= gap) val = 0;

			image[i * ndets + j] = val;
		}
	}

	for ( unsigned int i = 0; i < ndets; ++i ) {
		for ( unsigned int j = 0; j < ndets; ++j ) {
			fprintf(imafile, "%f ", image[i * ndets + j] * maxvalue);
		}
		fprintf(imafile, "\n");
	}
	fclose(imafile);

	for ( unsigned int i = 0; i < yres; ++i ) {
		for ( unsigned int j = 0; j < yres; ++j ) {
			int I = i * ndets / yres, J = j * ndets / yres;
			image1[i * yres + j] = image[ I * ndets + J ];
		}
	}
	glClear(GL_COLOR_BUFFER_BIT);		
	glRasterPos2f(-1, -1);
	glDrawPixels(yres, yres, GL_LUMINANCE, GL_FLOAT, image1); 

	sprintf(buffer,"%d", 0); Message(buffer, false);
	if (maxvalue > 1) {
		sprintf(buffer,"%d", (int)maxvalue); Message(buffer, true); 
	} else {
		sprintf(buffer,"%4.4f", maxvalue + 0.0005);
		for(int s = 0; s < 5; s++) buffer[s] = buffer[s+1];
		buffer[4] = '\0';
		Message(buffer, true);
	}
	glRasterPos2f(-1, -1);
	glReadPixels(0, 0, Xres, yres, GL_LUMINANCE, GL_FLOAT, image); 


	for ( int i = yres- 1; i >= 0; --i ) {
		for ( unsigned int j = 0; j < Xres; ++j ) {
			if (bar && j >= yres - 3 && j < yres) image[i * Xres + j] = (double)(i - 1)/yres;
			PrintTGAColor( tgafile, (double)image[i * Xres + j], 0.7 );
		}
	}
	fclose(tgafile);
}

bool compute = false;

void onDisplay( ) {
	glClearColor(0, 0, 0, 1);
	glMatrixMode(GL_MODELVIEW); glLoadIdentity();
	glMatrixMode(GL_PROJECTION); glLoadIdentity();

	if (!compute) {
		Compute();
		compute = true;
	}
}


//--------------------------------------------------------
// The entry point of the application 
//--------------------------------------------------------
int main(int argc, char *argv[]) {
	glutInit(&argc, argv);                // GLUT initialization
	glutInitWindowSize(256, 128);         // Initial resolution of the application Window is 600x600 pixels 
	glutInitWindowPosition(100, 100);     // Initial location of the application window
	glutInitDisplayMode(GLUT_RGB);        // Image = 8 bit R,G,B + double buffer + depth buffer
	glutCreateWindow("2dPET");        	  // Application Window is born
	glutDisplayFunc(onDisplay);           // Register event handlers 
	glutMainLoop();                    	  // Event loop  
	return 0;
}